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The Procellarum region is a broad area on the nearside of the Moon that is characterized 
by low elevations 1 , thin crust 2 , and high surface concentrations of the heat-producing 
elements uranium, thorium, and potassium 3 ' 4 . The Procellarum region has been interpreted 
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as an ancient impact basin approximately 3200 km in diameter ' , though supporting 
evidence at the surface would have been largely obscured as a result of the great antiquity 
and poor preservation of any diagnostic features. Here we use data from the Gravity 

o 

Recovery and Interior Laboratory (GRAIL) mission to examine the subsurface structure 
of Procellarum. The Bouguer gravity anomalies and gravity gradients reveal a pattern of 
narrow linear anomalies that border the Procellarum region and are interpreted to be the 
frozen remnants of lava-filled rifts and the underlying feeder dikes that served as the 
magma plumbing system for much of the nearside mare volcanism. The discontinuous 
surface structures that were earlier interpreted as remnants of an impact basin rim are 
shown in GRAIL data to be a part of this continuous set of quasi-rectangular border 
structures with angular intersections, contrary to the expected circular or elliptical shape 
of an impact basin 9 . The spatial pattern of magmatic-tectonic structures bounding 
Procellarum is consistent with their formation in response to thermal stresses produced by 
the differential cooling of the province relative to its surroundings, coupled with magmatic 
activity driven by the elevated heat flux in the region. 


The Procellarum KREEP Terrane (PKT) is defined by higher than average values of the 
surface abundances of potassium, rare earth elements, and phosphorus 3 ' 10 (Fig. 1). The PKT 
likely experienced a geodynamical history that differed from that of the rest of the Moon because 
of the elevated heat flow resulting from the high crustal concentrations of heat-producing 
elements 10 ’ 12 . The region encompasses the majority of the mare basalt provinces, including many 
that are not associated with known impact basins. The interpretation of the region as an impact 
basin was based on its distinctive composition and generally low elevation, together with the 
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photogeological interpretation of features as fragments of circular basin rings 5 ' 7 ’ 13 . The most 
prominent candidate ring structures are the mare shorelines and scarps on the western edge of 
Oceanus Procellarum and the northern edge of Mare Frigoris 5 (Fig. la; Extended Data Fig. 1). 
However, these arcuate segments span only a fraction of the circumference of the proposed 
basin, requiring that much of the proposed topographic rim was destroyed or modified beyond 
recognition. 

In this study, we use data from NASA’s GRAIL mission 8 to examine the subsurface 
structure of the Procellarum region. Bouguer gravity data (gravity field corrected for the 
contributions of surface topography) and gravity gradients (second horizontal derivatives of the 
Bouguer potential 14 ) reveal a distinctive pattern of anomalies surrounding the region (Fig. lb-c). 
These narrow belts of negative gravity gradients and positive gravity anomalies indicate narrow 
zones of positive density contrast in the subsurface. Previous work revealed a population of 
narrow, randomly oriented, ancient igneous intrusions that lack surface expressions 14 . In 
contrast, the PKT border anomalies are broader features that are spatially associated with the 
maria and appear to be part of an organized large-scale structure. These anomalies are the 
dominant features not associated with impact basins in the global gravity gradients, but only a 
portion of the western border anomalies in Oceanus Procellarum were noted in earlier gravity 
studies 15 . 

To investigate the source of the anomalies, we first inverted the gravity field in the 
spherical harmonic domain under the assumption that the anomalies arise from variations in the 
thickness of both the maria and the underlying feldspathic crust (see online-only Methods for 
details). We focus here on two models to illustrate the range of solutions: the first imposes an 
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isostatic condition on the pre-mare crust, and the second forces the amplitude of the relief along 
the mare-crust and crust-mantle interfaces to be equal and opposite in magnitude. For these two 
models, the average structure across two of the border anomalies at the northwest corner of the 
PKT suggest the presence of elongated mare-filled depressions in the feldspathic crust having 
widths of —150 km and depths of 2-4 km, and underlain by crust-mantle interfaces that are 
shallower than adjacent areas by 3-6 km (Fig. 2e-h; Extended Data Figs 2-3). If we instead 
assume that the PKT border anomalies arise from igneous intrusions in the subsurface 14 , 
inversions of the average gravity profiles across these two anomalies yield widths of 66 /g and 
82/jg km and vertical extents of 8/[ and 6), km for intrusions with elliptical cross-sections, 
assumed density contrasts of 550 kg/m 3 , and bottom depths of 25 km (Fig. 2c-d; see Methods). 

The spherical harmonic inversion solutions are consistent with thickening of the maria over 
linear depressions formed by crustal thinning, as could occur in volcanically flooded rift 
valleys 16 . The branching of anomalies of the western border structure and the triple-junction 
intersections at some corners are consistent with the attributes of planetary rifts. This 
interpretation is also supported by the broad elongated depressions surrounding the border 
anomalies beneath Mare Frigoris and western Mare Tranquilitatis, and the scarps found in the 
highlands adjacent to some of the border anomalies 5 . The inferred crustal thinning could arise 
from extension of the crust by 8-18 km (Extended Data Table 1). For the intrusion models, the 
large widths of the inferred intrusions (greatly exceeding the vertical dimensions), and the 
association of the gravity anomalies with mare basalts at the surface, suggest that dike-like 
intrusions are not solely responsible for the anomalies. A combination of crustal thinning, mare 
thickening, and intrusion by dike swarms provides the most likely explanation for the anomalies. 
The elevated heat flux in the PKT 10 coupled with passive mantle upwelling during rifting would 
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have led to widespread partial melting of the underlying mantle 16 , so extensional tectonics would 
have been accompanied by dike intrusion and volcanism. These dikes may represent the magma 
plumbing system that provided conduits connecting deep magma reservoirs to many of the 
nearside maria. 

The PKT border structures are the only known lunar structures consistent with large-scale 
rifting of the crust, a process more common on Earth, Venus, and Mars. The surface exposures 
of the maria overlying the border structures fonned 3.51±0.25 billion years ago (Ga; area- 
weighted mean and standard deviation) 17 , representing the final stages of the volcanic infilling of 
the structures. In contrast, the rest of the nearside maria exhibit a range of ages of 1. 2-4.0 Ga. 
Volcanic infilling of the rifts may have been a self-limiting process because the flexural response 
to the loading would have caused compression in the upper lithosphere, possibly closing off the 
magma conduits. This inference is supported by the observation of wrinkle ridges overlying and 
parallel to the border structures. Parallel wrinkle ridges flanking the Mare Frigoris border 
structure may also reflect structural control of the wrinkle ridges by buried tectonic structures. 

In a polar projection centred on the PKT, the border structures delineate a quasi- 
rectangular shape (Fig. 3). The arcuate scarps at the edges of Maria Frigoris and Procellarum 
that were previously interpreted as rim segments of a Procellarum basin are seen in the GRAIF 
data to be a small fraction of this continuous set of well-expressed structures that trace out a 
polygonal pattern consisting of predominantly straight sides and angular intersections (Extended 
Data Fig. 1). The northeast and northwest corners of the structure deviate from the proposed 
circular rim 5 by 215 km and 175 km, respectively. Only the discontinuous and poorly expressed 
anomalies in the southwest portion of the region are compatible with a circular rim. This quasi- 


112 

113 

114 

115 

116 

117 

118 

119 

120 

121 

122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 


rectangular pattern is in contrast with the circular or elliptical shapes of all other large impact 
basins 9 , including the hemisphere-scale Borealis basin on Mars, for which a continuous elliptical 
basin rim can be traced in topography and gravity data 18 . The interpretation of the PKT border 
structures as the rim of an impact basin would require hundreds of kilometres of horizontal 
deformation with large strain gradients to produce the angular corners. No evidence of such 
large-magnitude strain exists on the Moon 19 . Furthermore, the negative gravity gradients of the 
border structures do not match the signatures of known impact basins, such as the Imbrium and 
South Pole-Aitken basins, which are characterized by paired positive and negative gradients of 
equal amplitude flanking the rims and negative gradients throughout the basin interiors. 
Although it is not possible to disprove the existence of an ancient degraded Procellarum basin 
that lacks a clear geophysical signature, the geometry and gravitational signature of the structures 
bordering the PKT do not support the interpretation that they mark the rim of a basin. 

The formation and geometric pattern of the PKT border structures require an explanation. 
Although the gravity anomalies are consistent with either lava-flooded rift valleys or dense 
swarms of dikes, both interpretations require substantial extension across the border structures. 
The location of the structures at the edge of the PKT suggests that the elevated heat flux in this 
region 10 may have played a role in the extension inferred from the gravity modelling. In a state of 
thermal equilibrium, both the temperature and the rate of change in temperature in the 
lithosphere would be linearly proportional to the concentration of heat-producing elements in 
and/or beneath the crust. Thus, although the PKT was always warmer than its surroundings due 
to the high concentrations of heat-producing elements, it would have cooled at a greater rate due 
to declining radiogenic heat production 10 . The cooling lithosphere would then have experienced 
thermal contraction, which in turn would have caused horizontal extension at the margins. 
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Cooling by 600 K across a region 2000 km wide would have induced the equivalent of ~8 km of 
extension. We tested this hypothesis with a simple model of the thermal evolution and resultant 
stresses (see Methods). A finite difference model was used to represent the conductive thermal 
evolution of the Moon, given the equivalent of 10 km of KREEP basalt at the base of a 40-km- 
thick crust within a spherical cap 2000 km in diameter 10 ’ 11 . The model predicts a temperature 
decrease of the PKT relative to its surroundings of >600 K between 4.0 and 3.0 Ga, with the 
maximum cooling at the base of the crust (Fig. 4a; Extended Data Figs 4-5). 

The stresses resulting from the thermal contraction of the lithosphere between 4.0 and 3.0 
Ga were calculated with an elastic finite element model 20 . The far-field stresses on the opposite 
side of the planet were subtracted in order to isolate the effects of the PKT, since the mean stress 
in the lithosphere may have been affected by global contraction or expansion 14 ' 21 . Cooling and 
contraction of the lower lithosphere within the PKT caused extension, which induced 
compression in the elastically-coupled upper lithosphere inside the PKT, and extension 
throughout the lithosphere at the edge of the PKT (Fig. 4c). Similar results were obtained if the 
KREEP-rich material was distributed throughout the crust (Extended Data Figs 6-7). This 
extension may have been augmented by an early period of global expansion 14 . 

Many of the maria not associated with impact basins are found over the PKT border 
structures. Rise of magma to the surface in dikes requires the most tensile stress to be horizontal, 
as well as a vertical gradient in stress conducive to magma ascent 22 . The model predicts that the 
extensional zone bordering the PKT was conducive to magma ascent in dikes (Fig. 4d). In 
contrast, compressional stresses in the upper lithosphere within the centre of the PKT would tend 
to inhibit the rise of magma, except where this stress field was modified by later processes such 
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as impacts or loading and flexure of the lithosphere, or where magma ascent was aided by 
volatile exsolution or a pressurized magma chamber. 

In order to form the observed rectilinear pattern of structures, it is necessary to break the 
azimuthal symmetry assumed in the model. Volumetric contraction beneath a free surface 
generates fracture patterns with characteristic comer angles of 120°. This pattern results in six- 
sided polygons at scales ranging from 1-100’s of cm (e.g., mud cracks, columnar joints in 
basalt), to 1-100 m (e.g., thermal contraction polygons in permafrost), to 10 km (e.g., polygons 
from sediment compaction in the lowlands of Mars 23 ). However, as the size of the structure 
becomes large relative to the radius of the planet, surface curvature becomes important. A 
polygon with 120° corner angles will have five or four sides when the lengths of the sides reach 
32° or 80° of arc, respectively. The mean length of the PKT border structures is 2150 km or 71°, 
and the angles of the vertices range from 109° to 125°. Thus, at the scale of the PKT, a set of 
linear rifts intersecting at 120°-angle junctions around a contracting cap may result in a quasi- 
rectangular structure. 

We note a similarity in the pattern of structures to the south polar terrain (SPT) of Saturn’s 
icy moon Enceladus (Fig. 3; Extended Data Fig. 8) 24 ' 25 . Both the PKT and SPT are bordered by 
quasi-rectangular sets of tectonic belts with angular intersections that sometimes take the form of 
triple junctions. Both structures enclose regions approximately 70-80° in diameter of low 
topography 1 ' 25 , enhanced volcanic activity 10 ' 24 , and strongly elevated heat flow 10 ' 26 . However, we 
emphasize that there are important differences between the specific processes at work and the 
evolutionary histories of these two very different terrains - including the tidal source of the heat, 
the prevalence of compressional tectonism" ' , the likelihood of a subsurface ocean" , and the 
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possibility of a mobile lithosphere 28 on Enceladus. Nevertheless, the gross morphological and 
geophysical similarities between the PKT on the Moon and the SPT on Enceladus suggest the 
possibility of broad parallels in their geodynamic evolution, and that similar parallels may exist 
with other magmatic-tectonic centres (e.g., the northern lowlands of Mercury, an irregular 
depression -80° in diameter 29 that has experienced widespread volcanic resurfacing 10 ). 


METHODS SUMMARY 

Gravity gradients were calculated from the GRAIL Extended Mission gravity model 
GRGM900b with the SHTOOLS toolkit (available on-line at shtools.ipgp.fr). Spherical 
hannonic modelling of the mare and crustal thicknesses utilized SHTOOLS. The thermal models 
used a finite difference model of heat transfer for assumed concentrations of heat-producing 
elements in the PKT, crust, and mantle. Finite element modelling was conducted with the 
TEKTON code in an axisymmetric spherical geometry to examine the elastic stresses in the 
lithosphere. 


Online Content Any additional Methods and Extended Data display items are available in the online version of the 
paper; references unique to these sections appear only in the online paper. 
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Figure 1. Global maps of (a) topography, (b) Th concentration, (c) Bouguer gravity 
anomaly, and (d) gravity gradients on the Moon. All maps are simple cylindrical projections 
centred on the nearside. The circular rim of the proposed Procellarum impact basin 3 (black 
dashed line), the outline of the maria (white lines 17 ), and the extent of the PKT (red line, 
corresponding to a Th concentration of 3.5 ppm; ref. 4) are shown in a. Features discussed in the 
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text are labelled in a. 
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Figure 2. Gravity and subsurface structure of the PKT border structures, a-b, Maps of the 


modelled thickness of the maria and underlying feldspathic crust, c-d, Average Bouguer gravity 


profiles perpendicular to border anomalies 1 and 2 (see panel a for locations), e-h. Average 


cross-sections of the model results orthogonal to the border anomalies showing the mare (dark 


grey) and feldspathic crust (light grey) for two different sets of filters. The models in e-f impose 


the condition that the relief along the interfaces was in isostatic equilibrium prior to infilling by 


mare basalt, whereas the models in g-h impose the condition that the relief along the interfaces 
was equal and opposite in amplitude (see Methods for further details and Extended Data Fig. 3 


for results from additional models). 



232 Figure 3. Geometric pattern of the Procellarum KREEP terrane (PKT) border structures, 

233 with a comparison to the Enceladus south polar terrain (SPT). a, The border structures of 

234 the PKT highlighted by the gravity gradients trace out a quasi-rectangular pattern, enclosing b, a 

235 broad region of low elevations 1 , c, The SPT is similarly a region of low elevation 23 and high heat 

236 flow 26 (Extended Data Fig. 8) surrounded by a quasi-rectangular pattern of border structures. All 

237 maps are in a simple polar projection. In all panels, the circle corresponds to an angular diameter 

238 of 180° of surface arc, divided into 10° increments. 
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Figure 4. Predicted temperature and stress for the Procellarum region, a, Predicted 
temperature change of the PKT relative to its surroundings between 4.0 and 3.0 Ga. The 
Procellarum region is centred on the pole on the left side of the figure. The black line denotes 
the area expanded in panels c and d. b, In-plane horizontal elastic stress radial to the centre of 
the PKT at the surface predicted by the finite element model (where positive stresses are tensile; 
the far-field stress profile has been subtracted to calculate the relative stresses), c, Cross-section 
of the in-plane horizontal elastic relative stress, d, Predicted zones of magma ascent; dark grey 
indicates horizontal extension conducive to vertical dike formation, light grey indicates both 
horizontal extension and a vertical stress gradient more favourable to magma ascent than in the 
lithosphere far from the PKT, and red indicates areas in which magma will rise unassisted by 
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other factors. Cross-hatching indicates regions in which none of the criteria for magma ascent are 
met. The temperatures in a and stresses in b,c are both taken relative to the far-field values in the 
opposite hemisphere. 
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Methods 


Gravity gradients 

Gravity data was analyzed from gravity model GRGM900b obtained from observations 
during the primary and extended GRAIL missions. The Bouguer gravity anomaly model was 
generated for an assumed a crustal density of 2550 kg/m 3 (ref. 2). The Bouguer gravity gradients 
were calculated in the spherical harmonic domain 31 using the software archive SHTOOLS (freely 
available on-line at shtools.ipgp.fr). The eigenvalues of the horizontal gravity gradient tensor 
(Gn, G 22 ), representing the values of the maximum and minimum curvature of the potential field 
at each point, were then calculated. As was done previously 14 , the eigenvalues were combined 
into a single value (the maximum-amplitude horizontal gradient, or G /,/,) representing the second 
horizontal derivative of maximum amplitude at each point on the surface: 


r n if r n 

> 

r 

a 22 1 

r 22 if r n 

< 

F 

a 22| 


where |x| indicates the absolute value of x. This maximum-amplitude horizontal gradient 
represents the gradient orthogonal to any structures that dominate the local gravity, regardless of 
their orientation. The gravity gradients are given in units of Eotvos (1 E = 10" 9 s' 2 ). The gravity 
gradients were used to reveal the presence of discrete subsurface structures, whereas the Bouguer 
gravity anomaly and potential were used in all subsequent analyses. 

In this representation of the gravity gradients, a positive density anomaly will produce a 
negative gravity gradient, whereas a step function density anomaly will produce a symmetric pair 
of positive and negative gravity gradients flanking the step. For this reason, the mantle uplift 
beneath large impact basins is expressed as an outer ring of positive gravity gradients and an 
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inner ring of negative gravity gradients. Thus, although some of the border structures are near 
the edges of the overlying maria, the gravity gradient signatures are not consistent with the 
anomalies expected to arise from edge effects of the maria. Furthermore, the northern border 
anomaly is approximately centred within Mare Frigoris, and the western border structure exhibits 
three branches that are offset from the edge of the overlying Oceanus Procellarum by as much as 
600 km. The average Bouguer gravity profiles perpendicular to the border structures reveal 
narrow positive Bouguer anomalies (Fig. 2c, d). The elongated negative gravity gradients and 
positive Bouguer gravity anomalies bordering the Procellarum KREEP Terrane (PKT) are most 
simply explained by elongated positive density anomalies. 

In previous work focusing on narrower structures in the lunar gravity gradient field 
interpreted as giant dikes or swarms of dikes, we calculated the gradients using a high-pass filter 
at degree and order 50, which emphasized shorter-wavelength structures 14 . The focus of the 
present work is on the longer-wavelength border anomalies surrounding the PKT, which have 
significant power at degrees less than 50. Thus, the gravity gradients were calculated between 
degrees 2 and 400, with a cosine-shaped taper applied between degrees 350 and 400. Two of the 
border anomalies in the northwest part of the region coincide with ancient igneous intrusions 
identified in the previous study of the short-wavelength gravity gradients 14 . However, the 
majority of the giant dikes identified in that study are narrower structures that lack a surface 
expression and are distributed randomly across the planet 14 . In contrast, the PKT border 
anomalies are longer-wavelength structures that occur within the maria and appear to be part of a 
large-scale organized structure. 

In order to highlight the true shape of the PKT border anomalies, the Bouguer gravity 
data and gradients were plotted in a simple polar projection, preserving the distance between 
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each point and the origin, and thus preserving the shape of features centred on the origin. The 
global Bouguer gravity gradient map in cylindrical projection (Fig. 1) appears to show a 
pentagonal structure encompassing the PKT. However, re-projection in a polar projection 
centred on the region (Fig. 3a) reveals that the structure as a whole is dominantly quasi- 
rectangular. The pentagonal appearance in the cylindrical projection is a result of both the 
distortions at high latitude in that projection and a kink in the northern border structure at its 
mid-point. 

A previous study 5 mapped possible ring structures associated with the Procellarum basin 
on a Lambert azimuthal equal-area map of the nearside of the Moon. A comparison of the 
GRAIL gravity gradients with this map (Extended Data Fig. 1) reveals that the majority of the 
mare shorelines and major scarps identified in that study parallel the Procellarum border 
anomalies, and a significant fraction of the wrinkle ridges overlie the border anomalies. 
However, the angular corners apparent in the gravity gradients are missing or rounded off in the 
mapped surface structures. The scarps and mare shorelines adjacent to the border anomalies are 
consistent with their interpretation as lava-flooded rifts, and the alignment of wrinkle ridges over 
the border anomalies is consistent with the flexural stresses expected to arise from the narrow 
loads inferred from the gravity data. The tracing of these structures on a Lambert azimuthal 
equal-area map, which does not preserve angles and causes significant distortions around the 
edges due to the non-linear radial distance scale, contributes to the apparent circularity of the 
structures. This distortion is particularly prominent for the northwest corner of the PKT border 
structures, which occurs near the limb of the Moon where the distortion is at its greatest. 
Nevertheless, even in this projection the border anomalies clearly delineate a polygonal structure. 
A simple polar projection centred on the Procellarum region preserves the distance from the 
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center to all points and thus provides a more accurate depiction of shapes centered on the origin. 
Only the discontinuous structures in the southwest corner of the Procellarum region are 
consistent with a circular pattern. 

Gravity inversions 

Long-wavelength Bouguer gravity anomalies on the Moon are thought to arise largely 
from variations in the relief along the crust-mantle interface 2 ' 32 . In contrast, because the 
gravitational potential of short-wavelength anomalies attenuates rapidly with elevation, most of 
the observed high-degree power in the Bouguer gravity must arise at depths shallower than the 
crust-mantle interface. At intermediate degrees, the origin of the gravity anomalies depends on 
the geodynamic setting. For the case of the PKT, the vast majority of the border anomalies occur 
beneath maria, and thus the anomalies likely arise at least in part from variations in the relief 
along the mare-crust interface. However, some minor branches extend off from the main border 
anomalies into the surrounding crust outside the maria, suggesting that at least some component 
of intrusive dikes and/or uplifted crust-mantle interface contributes to the anomalies. We 
consider both possibilities in our analysis. 

The width of the gravity anomalies and their association with mare basalts at the surface 
suggest that the anomalies may be the result of local thickening of the maria above linear 
tectonic structures and/or uplift of the crust-mantle interface beneath those structures. To 
investigate this scenario, we inverted the gravity data in the spherical harmonic domain by 
downward continuing the Bouguer gravity to the appropriate radii and iteratively solving for the 
spherical harmonic coefficients describing the relief along the density interfaces of interest, 
taking into account the finite-amplitude effects of that relief 52 . This approach has been applied 
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previously for calculating the relief along the crust-mantle interface 2 ’ 2 , but here we wish to solve 
for the relief along both the mare-crust and crust-mantle interfaces. We first calculated the 
Bouguer gravity using the density of mare basalt, since the maria comprise the top layer in our 
three-layer model (mare, crust, and mantle). We adopt a mare density of p,„=3150 kg/m 3 , based 
on the average of measured densities of Apollo mare samples 3 ’. The Bouguer anomaly was then 
used to calculate the relief along the mare-crust and crust-mantle interfaces. 

The solution for the relief along two different subsurface density interfaces is inherently 
non-unique. In order to capture a range of possible solutions, we consider different filters to 
parse the gravity anomalies between the crust-mantle interface and the mare-crust interface. We 
designed a filter wi to allow us to specify the desired ratio, /, between the relief along the crust- 
mantle interface and that along the mare-crust interface, taking into account the degree- 
dependent amplification of the gravity anomalies during their downward continuation to the 
mean depth of the interface of interest: 


{RJR„TXp„-p,)+{ r m !R«T(pm-p} f 

where l is the spherical harmonic degree, p c is the density of the feldspathic crust, p M is the 
density of the mantle, p m is the density of the mare, Ro is the mean planetary radius (1737.15 km; 
ref. 1), R m is the mean radius of the mare-crust interface, and Rm is the mean radius of the crust- 
mantle interface. This filter was applied in calculating the relief along the crust-mantle interface, 
and the remaining Bouguer gravity was then used to calculate the relief along the mare-crust 
interface. We assumed densities of 2550 kg/m 3 and 3220 kg/m 3 for the crust and mantle, 
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respectively, on the basis of previous GRAIL analyses 2 . We assumed a mean radius of the crust- 
mantle interface of 1697.15 km, resulting in a mean crustal thickness of 40 km, and a mean 
radius of the mare-crust interface of 1736.15 km. The filters used for the models depicted in Fig. 
2 are shown in Extended Data Fig. 2. The first model represents the case in which the feldspathic 
crust was in a state of isostasy prior to infilling by the mare, leading to a ratio / of pJ(pM-Pc ) 
(Extended Data Fig. 2a). In this model, isostasy is defined using the simple criterion of equal 
masses in adjacent columns. If some of the volcanic infilling of the structures occurred in 
parallel with the extensional tectonics, the resulting load would have driven added subsidence, 
which would have changed the ratio between the relief along the mare-crust and crust-mantle 
interfaces. The second model represents the case in which the relief along the two interfaces was 
equal and opposite in amplitude, with/taking on a value of 1 for degrees >10. However, because 
the long-wavelength topography of the Moon is largely isostatic, we assumed the isostatic ratio 
for / for degrees 1-3, with a linear transition between the isostatic and equal amplitude values 
over degrees 3-10, and the equal amplitude value from degrees 10 to 125 (Extended Data Fig. 
2b). These two models serve to illustrate the range of possible solutions and the relative 
insensitivity of the inferred extension to the model assumptions. A low-pass cosine taper from 
degrees 125 to 150 was applied to all models. 

The resulting models match the gravity data but do not take into account the effects of 
flexure, which would perturb the interface depths relative to their elevations prior to mare 
loading and thus alter the assumed pre-loading ratio between the relief along the interfaces. 
Although the models were applied globally, the results are not valid in areas outside the maria. 
Similarly, crustal thickness models that neglect the high density of the mare basalt and the 
possible variations in mare thickness will have errors within the maria. The mean radius of the 
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mare-crust interface was chosen so as to bring the base of the maria within the Procellarum 
region below the surface over most of the observed maria. However, the modeled long- 
wavelength variations in the thickness of the maria are poorly constrained because of the 
ambiguity between the gravitational effects of variations in the relief along the mare-crust and 
crust-mantle interfaces. As a result, the distribution of areas with predicted mare thicknesses 
greater than zero only approximately matches the observed distribution of the maria. 
Nevertheless, the shortwavelength variations in the thickness of the maria beneath the border 
anomalies are robust, given the model assumptions. 

The density of the lunar mantle beneath the PKT is not known. The process responsible 
for concentrating the KREEP-rich materials on the nearside of the Moon may have also brought 
dense ilmenite-rich cumulates to the base of the crust on the nearside’ 4 . Overturn of the 
buoyantly unstable magma ocean cumulates would have mixed this material to deeper levels in 
the lunar mantle’ 3 ' 36 , but this overturn is limited by the high viscosity of the solid ilmenite-rich 
cumulates and will only occur for a limited range of scenarios 37 . It is possible that a mixture of 
olivine and ilmenite-rich cumulates sa nk as solid diapirs, leaving behind a portion of the 
ilmenite-rich material at shallower levels 37 . To account for the possibility of shallow ilmenite- 
rich material beneath the PKT, we also consider a high-mantle-density endmember model with 
an assumed mantle density of 3500 kg/m 3 , representative of the density of the late stage 
crystallization products from the magma ocean’ 6 . The higher mantle density reduces the 
predicted mantle uplift beneath the border structures, and similarly reduces the predicted 
extension. 

We also considered two additional end-member scenarios in our gravity models. For one 
model, we assumed that all of the gravity anomalies at degrees >10 arise from variations in the 
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thickness of the maria. This model required a mean mare-crust interface radius of Ro-6 km in 
order to bring the mare-crust interface below the surface in the regions of interest. For another 
model, we assumed that all of the gravity anomalies at degrees >10 arise from variations in the 
relief along the crust-mantle interface. This model became unstable at higher degrees due to the 
amplification of the high-degree gravity anomalies during downward continuation to the mean 
depth of the crust-mantle interface, so a cosine taper was applied between degrees 75 and 100 to 
stabilize the solution. As a result, this model is a factor of 1 .6 coarser in resolution than the other 
models. This result provides further evidence that the short-wavelength gravity anomalies must 
arise from density anomalies at depths more shallow than the crust-mantle boundary. This model 
ascribing all of the Bouguer gravity to variations along the crust-mantle interface is comparable 
in resolution to the global GRAIL crustal thickness models 2 (low-pass filtered with an amplitude 
of 0.5 at degrees 87 and 80, respectively, corresponding to spatial wavelengths of 63 and 68 km). 
In contrast, the models ascribing a substantial fraction of the Bouguer gravity to the shallower 
mare-crust interface are higher in resolution (low-pass filtered with an amplitude of 0.5 at degree 
137, corresponding to a spatial wavelength of 40 km). For both models described above, we 
assumed that variations in the the top and bottom surfaces of the feldspathic crust from degrees 1 
to 3 were isostatically compensated before mare flooding, with a linear transition to the desired 
filter from degrees 3 to 10. These final two models are not likely to be accurate representations 
of the subsurface structure, but they bracket the range of possible solutions. 

The predicted relief along the interfaces was used to calculate the thic kn esses of the crust 
and maria (Extended Data Fig. 3). The broad patterns of mare thickness in this region as 
indicated by the models are highly uncertain, due to the non-uniqueness of the division of the 
gravity anomalies between the mare-crust interface and the crust-mantle interface. In some 
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areas, the predicted base of the mare rises above the surface, indicating the need for subsurface 
mass deficits as could arise from additional variations in the crustal thickness or density in order 
to explain the observed gravity within the context of this model. These errors outside the maria 
do not affect the predictions for the Procellarum border structures. The local thickening of the 
mare over the western Procellarum border structure is to first order consistent with maps of the 
mare thickness derived from geological constraints, such as the burial depths of impact craters, 
which show local thickenings of up to >1.5 km along this structure 38 . Models combining the 
effects of dikes with the relief along the mare-crust and crust-mantle interfaces would predict 
narrower dikes than models that ascribe the entire gravity anomaly to the presence of dikes, and 
reduced relief along the density interfaces relative to models without dikes. 

The extension across the structures was calculated from the thickness of the feldspathic 
crust by integrating the fractional crustal thickness anomaly across the structures: 

where A L is the change in length, c(x) is the thickness of the feldspathic crust as a function of 
location, and Co is the mean thickness of the crust on either side of the structure. The extension 
was calculated between the shoulders on either side of the rift for each model, encompassing a 
zone 131 and 152 km wide for anomalies 1 and 2, respectively. The calculated extension and 
corresponding extensional strain across the structures for each of the models are given in 
Extended Data Table 1. The models with an isostatic ratio between the relief at the top and 
bottom of the feldspathic crust predict greater extension because a larger fraction of the gravity 
signal is downward continued to the crust-mantle interface, resulting in greater amplification of 
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the short-wavelength anomalies. The extension calculated using the crustal thickness models is 
an upper bound because some contribution to the gravity anomaly arising from the mechanical or 
thermal reduction of the crustal porosity beneath the mare load and surrounding the intruded 
dikes is likely. 

We next inverted the Bouguer gravity over the PKT border structures for the best- fit 
dikes using a Monte Carlo approach. The sources of the anomalies were represented as density 
anomalies with elliptical cross-sections in the vertical plane perpendicular to the long axes of the 
anomalies, of assumed density contrast and bottom depth and unknown width and top depth. The 
bottom depths were set to the typical crustal thickness within the PKT of ~25 km (ref. 2), and the 
density contrasts were set to 550 kg/m 3 , corresponding to a crustal density of 2550 kg/m 3 (ref. 2) 
and an intrusion density of 3100 kg/nT (ref. 33). Dikes with elliptical cross-sections were then 
constructed from a large number of rectangular prismatic elements, and the gravity anomaly was 
calculated from those prisms 39 . The best-fit solutions were found using a simple Markov chain 
Monte Carlo (MCMC) approach 14 . The one-standard-deviation (1-a) confidence intervals on the 
best-fit solutions were obtained by using a Metroplis-Hastings MCMC to test 20,000 models and 
analyzing the histograms of the resultant model parameters 14 . If the volume of the dike is 
accommodated solely by horizontal extension, then the resulting extensions for anomalies 1 and 
2 are 21 km and 20 km, respectively, given intrusion into a 25-km-thick crust. 

Thermal modelling 

The thermal evolution of the PKT was modelled following earlier work by Wieczorek 
and Phillips 10 and Grimm 11 , under the assumption of conductive heat transfer through the 
mantle. The results of this work are primarily sensitive to the temperatures in the lithosphere, 
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which are dominated by the concentration of heat-producing elements in the crust and the 
conductive heat transfer through the lithosphere. Although early convection beneath the PKT 
was possible 12 , this convection would have had only a second-order effect on the temperatures in 
the lithosphere. We used a finite difference approach to solve the spherical axisymmetric thermal 
diffusion equation. The model was benchmarked against the analytic solution for half-space 
cooling from an instantaneous temperature change applied to the surface, as well as by 
comparison with the results of previous work 10 . The model cells were divided into crust, mantle, 
and KREEP components. 

The PKT was represented by a spherical cap 2000 km (66°) in diameter in which the 
concentration of heat-producing elements was enhanced. The lack of similarly high 
concentrations of heat-producing elements on the farside is supported by the lack of evidence for 
KREEP-rich material within or surrounding the South Pole-Aitken impact basin 3 . The cause for 
this concentration of incompatible elements on the nearside is not known, but it may be related to 
a degree- 1 Rayleigh-Taylor instability that arose from the gravitational instability of the dense 
ilmenite rich cumulates formed in the late stages of magma ocean crystallization’ 4 . The crustal 
thickness was set to a unifonn value of 40 km in order to isolate the effect of the concentration of 
heat-producing elements in the PKT. The effect of the thicker crust outside of the PKT is less 
than the uncertainties in the concentration of heat-producing elements and the thermal 
conductivity of the crust and PKT. We assumed a thermal conductivity of 2.0 W/m K for the 
crust and KREEP-rich material, and 3.0 W/m K for the mantle. The densities of the crust/PKT 
and mantle were set to 2550 and 3200 kg/m 3 , respectively, and a specific heat of 1200 J/kg K 


was assumed for all materials. 
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Previous studies favoured a 10-km-thick layer of KREEP basalt at the base of the 
crust 10 ' 12 , but other workers have argued that this scenario is not compatible with the gravity and 
topography of the region and generates too much melt 11 . In our nominal model, we included a 
10-km-thick layer of KREEP basalt at the base of the crust. We also considered the case of a 10- 
km-thick layer of KREEP basalt distributed uniformly throughout a 40-km-thick crust. We 
assumed a U concentration in the KREEP basalt of 3.4 ppm by weight, and concentrations of 
0.14 ppm and 6.8 ppb in the crust and mantle, respectively 10,12 . We assumed a K/U ratio of 2500 
and Th/U ratio of 3.7 in all materials 12 . The enhanced concentration of KREEP is given an abrupt 
edge in the thermal model for simplicity. The thermal effects of this edge are broadened over the 
thermal diffusion length scale (-50 km for 100 million years), whereas the stress effects are 
spread out over a distance comparable to the flexural half-wavelength (-540 km for a lithosphere 
thickness of 50 km). The overall stress pattern would be unaffected by tapering the margins of 
the KREEP terrane over length scales of this order. The effects of melting and melt extraction 
on the temperature evolution were neglected. Extraction of melt would reduce the magnitude of 
the thermal anomaly in early time steps and decrease the amount of cooling by a modest amount, 
but would not change the character of the results. 

High temperatures throughout the lunar interior are expected after accretion and 
solidification of the magma ocean 36 . The model was initialized with an approximation of an 
adiabatic temperature gradient throughout the model domain 10 , increasing linearly from 1450 K 
at the surface to 1500 K at the core-mantle boundary at a radius of 438 km. This temperature 
profile represents the temperature at the end of an early convective period. In the absence of an 
early period of convection, the temperatures at the top of the mantle after magma ocean overturn 
would have been similar 36 . The top boundary condition was set to a constant temperature of 250 
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K, approximating the radiative equilibrium temperature of the lunar surface. A constant heat 
flux of 0 was applied as the basal boundary condition at the core-mantle boundary. The model 
begins at t=0 (4.5 Ga) and was run forward in time for 4.5 billion years. The change in 
temperature with time was calculated between 4.0 Ga (somewhat before the onset of the 
geological record) and 3.0 Ga, bracketing the period during which the majority of the maria 
formed 17 ' 40 ’ 41 . It is only the change in temperature that generates thermal stresses in the 
lithosphere, so even though the PKT was always warmer than its surroundings, its time evolution 
was characterized by net cooling and thermal contraction because it cooled at a faster rate. The 
temperature change of the PKT relative to the surroundings was also calculated for illustration 
purposes by subtracting the temperature profile at the antipode of the PKT. The absolute change 
in temperature was used in all stress modelling, but the relative temperature change serves to 
highlight the evolving thermal anomaly beneath the PKT. 

The changes in temperature as functions of time at 25-km-depth (the midplane of the 50 
kin-thick lithosphere assumed for the stress modeling) both within and outside of the PKT are 
shown in Extended Data Fig. 4. Both scenarios for the distribution of KREEP-rich material 
show similar patterns, but the model with an isolated KREEP-rich layer beneath the crust 
experiences an early phase of warming in the first few hundred million years. Between 4.0 and 
3.0 Ga, both models predict substantially more cooling in the PKT than elsewhere. The mantle 
immediately below the PKT follows a similar pattern of cooling with time as a result of the 
decline in heat production within the PKT. In contrast, the mantle at deeper levels warms up as 
it slowly comes into thermal equilibrium with the overlying KREEP material 10 . The changes in 
temperature of the upper and lower mantle with time are most pronounced for the case of 
KREEP at the base of the crust (Extended Data Fig. 5). However, the net effect of the cooling 
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upper mantle and warming lower mantle approximately cancel out. The temperature changes 
predicted here are somewhat larger than those of Wieczorek and Phillips 10 as a result of the 
different ratios between the concentrations of heat -producing elements 12 and the neglect of latent 
heat and melt extraction effects in this study. Reducing the concentration of radiogenic isotopes 
or taking into account melt extraction would reduce the magnitudes of the predicted temperature 
changes and stresses, but would not affect their spatial patterns. 

There is substantial uncertainty in the early thermal state of the Moon. The variation of 
temperature with depth after accretion and solidification of the magma ocean depends strongly 
on the timescale of accretion , the depth of the magma ocean - ’ and the possible gravitational 
overturn of the magma ocean cumulates 36 ' 37 . However, our models depend primarily on the 
temperatures within the lithosphere, which are dominated by the time evolution of the heat 
production within the crust. By 4.0 Ga, at which time we begin tracking the temperature changes 
to calculate the strain, the effect of the assumed initial condition on the temperatures in the 
lithosphere is greatly reduced. The early period of thermal equilibration of the lithosphere is 
reflected in the -200 million year period of increasing temperature for the case of KREEP-rich 
material concentrated at the base of the crust (Extended Data Fig. 4). The possible persistence of 
mantle convection throughout the time period of interest 12 would affect the distribution of 
temperature with depth in the mantle but would have little effect on the time evolution of the 
temperature in the lithosphere. 

Both Apollo seismic observations 43 and GRAIL gravity observations 2 indicate that the 
Moon’s upper crust is fractured and porous, possibly to a depth of -20 km. This porosity is 
likely to reduce the thermal conductivity of the upper crust 44 . The viscous closure of porosity is a 
thermally activated process 2 ' 45 , so the higher temperatures within the PKT (Extended Data Fig. 
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4) may have decreased the crustal porosity and increased the thermal conductivity in the PKT 
relative to its surroundings. This increased thermal conductivity would have acted to accelerate 
the cooling of the PKT relative to that shown in Extended Data Figs 4-5. We have not attempted 
to model this process in detail, but we note that it will positively reinforce the thermal evolution 
shown here. 

Stress modelling 

The stresses resulting from the changes in temperature with time were modelled using the 
Tekton finite element software 20 in a spherical axisymmetric geometry subject to a uniform radial 
gravitational acceleration. In order to provide adequate spatial resolution in the PKT, the model 
domain was limited to the elastic lithosphere, assumed to be 50 km thick (see discussion below). 
The bottom boundary condition represented the restoring force of the mantle with a pressure that 
varied with depth, whereas elements were free to move in both vertical and horizontal directions. 
The effects of the buoyant upward pressure arising from thermal anomalies in the mantle below 
the PKT were applied to the bottom boundary as an additional pressure term that varied with 
location on the basis of the thermal model results. This pressure term was calculated as the depth 
integral of the density contrast relative to background density, scaled by the gravitational 
acceleration. Although considerable thermal anomalies are predicted in the sub-lithospheric 
mantle beneath the PKT, the effects of the cooling upper mantle and warming lower mantle 
largely cancel out. The remaining pressure contributes to a broad upwarping of the surface 11 but 
has little effect on the short-wavelength stresses that are the focus of this analysis. The final 
topography and gravity anomalies over the PKT as a whole would have been strongly affected 
by the flexural resistance of the lithosphere 11 , the thinning of the crust within the PKT 2 , and 
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loading by the maria 12 . The excess basal pressure far from the PKT, representing the effects of 
net global expansion or contraction, was subtracted from the basal pressure condition throughout 
the model. The net volume change of the interior could add a uniform compressional or 
extensional horizontal stress to the lithosphere, depending on the early thermal history of the 
Moon 14 ' 21,42 . The model domain of a 50-km-thick lithosphere stretching from pole to pole was 
divided into 600 nodes in the azimuthal direction and 20 nodes in the radial direction, resulting 
in element dimensions of 9. 1 by 2.5 km at the surface. 

The predicted temperature change between 4.0 and 3.0 Ga (Extended Data Figs 5b, d) 
was used to calculate the resulting instantaneous elastic stresses in the model elements prior to 
any deformation 46 : 


a = a v AT 


E 

3(1 -2v) 


where o is the stress (taken here to be isotropic), a v is the volumetric coefficient of thermal 
expansion (assumed to be 2x10° K' 1 ), E is Young’s modulus (assumed to be 100 GPa, which is 
likely appropriate for the lower crust in which the greatest contraction occurs), and v is Poisson’s 
ratio (assumed to be 0.25). These pre-strain thermal stresses were added to the lithostatic stresses 
for the initial condition for the finite element model. Imposing the effects of thermal contraction 
with the pre-strain stresses allows the resultant deformation and its effects on the stress field to 
arise self-consistently in the model. 

The elastic stresses were calculated relative to the far-field values at the opposite side of 
the planet in order to isolate the effects of thermal contraction of the PKT. Geological and 
geophysical evidence suggests that the net stress state of the Moon may have evolved from 


692 

693 

694 

695 

696 

697 

698 

699 

700 

701 

702 

703 

704 

705 

706 

707 

708 

709 

710 

711 

712 

713 


global expansion and extension to contraction and compression over the course of its thermal 
evolution 14 ' 21 . In this scenario, the net global stress change at the time of fonnation of the border 
anomalies may have been small. However, theoretical models have shown that an early period of 
global expansion is difficult to generate for many likely lunar formation scenarios 42 . We put this 
question aside and focused instead on the local stresses within and surrounding the PKT relative 
to the typical stresses far from the region. These stresses would have been modified by the 
global stress state at the time of interest by the addition of a uniform compressional or 
extensional horizontal stress. In addition to the relative stresses, we also show the difference 
between the in-plane horizontal and vertical stresses (Oh-o v ) and the deviatoric horizontal stress 
(Oh-a p ), where a p is the pressure or mean stress value over all three directions (Extended Data 
Fig. 6). The width of the zone of predicted extension (-400 km) is somewhat wider than the 
observed border structures (-200 km), but localization of the strain release would likely have 
occurred if the structures are analogous to lava-flooded rifts. Similar stress patterns are predicted 
if KREEP-rich material is distributed uniformly through the crust, though the magnitudes of the 
stresses are reduced (Extended Data Fig. 7) because of the reduced temperature changes 
(Extended Data Fig. 5c, d). 

The stresses predicted by the model are dominated by the simple horizontal contractional 
stresses within the lithosphere. However, volumetric contraction does induce small changes to 
the surface topography, which generates bending stresses of small magnitude. Models in which 
the vertical displacement was set to zero at either the top or the bottom of the model domain 
resulted in similar stress fields, demonstrating that bending stresses do not contribute 
significantly. 
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The modelling in this study is intentionally simple in order to isolate the effect of the 
contracting cap within the PKT. This analysis did not consider the effects of spatial or temporal 
variations in the lithosphere thickness. If the base of the lithosphere follows the 800°C isotherm 
in the dry lunar mantle 47 , the lithosphere thickness would range from 33 km within the PKT to 
216 km outside of the PKT at 3 Ga for the case of KREEP concentrated at the base of the crust, 
with a smaller amplitude change in lithosphere thickness for KREEP distributed throughout the 
crust. The adopted lithosphere thickness value of 50 km is intermediate between the predicted 
values within and outside of the PKT. This assumption is appropriate, since the tensile stresses 
of interest are at the edge of the PKT. Since the dominant source of stress is the horizontal 
contraction of the lithosphere within the PKT, the stresses for the case of a variable lithosphere 
should be similar. Models that were run with a spatially variable lithosphere thickness predicted 
similar results but were subject to artefacts resulting from the large distortions in the element 
grid at the edge of the PKT. This model represented only the elastic stresses within the 
lithosphere. A viscoelastic model of the lithosphere and underlying mantle would predict a 
viscous transition zone at the base of the lithosphere within which the stresses decreased to zero 
at depth. Coupling of the thermal and viscoelastic evolution would result in a lithosphere that 
thickens with time, and would likely reduce the magnitude of the predicted extension, but would 
not change the character of the results. Within the PKT, the stresses are characterized by 
compression in the upper lithosphere and extension in the lower lithosphere, whereas at the 
edges of the PKT the extensional stress reaches the surface. However, the frictional strength at 
the surface should approach 0 MPa, allowing release of the shallow compressional stresses. 
Brittle compressional failure of the frictionally weak upper lithosphere throughout the PKT 
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would allow further contraction of the spherical cap, significantly enhancing the extension at its 
margins. 

In order to model directly the formation of the observed border structures, it would be 
necessary to localize the extension through tectonic failure. The localization of the extensional 
failure at discrete rift zones in the border structures would be dependent on the strain rate, 
rheology, and crustal thic kn ess 48 . Failure at the edges of the PKT would relieve the stresses in 
the interior and allow the spherical cap to pull away from the surrounding lithosphere. Future 
work is needed to model more directly the formation of these border structures. In this work, we 
simply show that thermal contraction of the PKT predicts extension at its edges, providing a 
straightforward model for generating the PKT border structures. Additional stresses arising from 
uplift or subsidence of the lithosphere 11 ’ 12 and magmatic processes would have also likely played 
a role. 

Zones favourable to the ascent of magma-filled dikes through the lithosphere were 
identified as those experiencing in-plane horizontal extension relative to the vertical stress and a 
favourable vertical stress gradient. Horizontal extension is required for the formation of vertical 
dikes, which would otherwise flatten out to produce horizontal sills. In addition, the upward 
propagation of the dikes requires that the vertical gradient in the confining horizontal stress in 
the lithosphere (dah/dz, where positive stresses are tensile, z is positive upward, and Oh includes 
both the lithostatic stress and the added tectonic stress) be greater than the hydrostatic pressure 
gradient in the magma, causing the lower tip of the dike to pinch shut while the upper tip 
propagates upward. The low density of the lunar crust 2 is an impediment to the rise of magma, 
even in a neutral stress state. Magma ascent is favoured in cases in which the upper lithosphere 
is in a state of extension relative to the lower lithosphere 49 . For a magma density of 2900 kg/m 3 , 
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this state corresponds to a vertical gradient in the horizontal stress in the lithosphere in excess of 
4.7 MPa/km. However, the stress gradient in the upper portions of a conductively cooling 
lithosphere with internal heat production and basal heating is generally not conducive to magma 
ascent as a result of the increasing horizontal extension with depth caused by the declining 
geothermal gradient in the lithosphere with time. This problem could be ameliorated by a failure- 
induced reduction in the extensional stresses in the lower crust, by volatile exsolution within the 
magma to enhance the driving force for magma ascent 50 , or by a pressurized magma reservoir at 
depth. We use the criterion of a vertical stress gradient >4.7 MPa/km for unassisted magma rise, 
and we also look at the stress gradient relative to the far-field value antipodal to the PKT to 
assess the relative tendency for magma to rise through the lithosphere if assisted by other factors 
as discussed above. 

By these criteria, the zone at the margin of the PKT experiences stresses most conducive 
to magma ascent and eruption. Extensional horizontal stresses radial to the centre of the PKT 
would facilitate the formation of circumferential dikes throughout the full vertical extent of the 
lithosphere. The stress gradient in this zone is more conducive to magma ascent than anywhere 
else on the Moon. For the case of heating by a layer of KREEP at the base of the crust, magma 
would be predicted to rise unassisted to the middle of the lithosphere, whereas further ascent 
would require additional driving forces (Extended Data Fig. 6c). For the case of heating by 
KREEP distributed throughout the crust, the zone at the edge of the PKT is still the preferred 
location of magma ascent, but some added driving force such as volatile exsolution or a 
pressurized magma chamber is required for the rise of magma into the crust (Extended Data Fig. 


7c). 
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The model also predicts changes to the surface topography. The thermal contraction of 
the lithosphere with time causes surface subsidence due to the vertical component of that 
contraction. Additionally, the horizontal shortening of the spherical cap centred on the PKT 
results in further subsidence due to the effects of the membrane stresses. Taking into account 
both the thermal contraction of the lithospheric cap in the PKT and the effects of thermal 
anomalies in the mantle, our models predict changes in surface topography less than 0.5 km 
during the period between 4.0 and 3.0 Ga. This result cannot be directly compared with the 
observed topography because it represents only the change in topography over a fraction of lunar 
history. However, we note that the predicted elevation changes are smaller than the observed 
relief. The topographic depression within Procellarum cannot be explained by thermal 
subsidence alone. 

Previous work found that the patterns of uplift predicted by thermal models of the PKT 
are difficult to reconcile with the observed long-wavelength gravity and topography 11 ’ 12 . 
However, the gravity and topography within the PKT are also likely affected by variations in the 
thickness of the crust and by possible density anomalies in the underlying mantle 12 . The low 
topography within the PKT may also be affected by a reduction in the porosity of the crust from 
thermal annealing 45 . For a lunar crustal porosity of 12% (ref. 2), annealing the pore space in the 
lower 10-20 km would reduce the surface elevation by 1.2-2. 4 km, consistent with the observed 
relief. 

We suggest that the quasi-rectangular pattern of border structures surounding the PKT is 
consistent with the intersection of linear rifts at 120°-angle triple junctions when the effect of the 
curvature of the surface is taken into account. Although the PKT border structures display some 
intermediate kinks and intersections, the overall planform is quasi-rectangular. Similarly, 
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although small-scale contraction crack polygons of all types in nature often have highly irregular 
forms, in the absence of competing effects (such as progressive subdivision of the polygons) the 
average structure is hexagonal because of the dominance of 120°-angle triple junctions at the 
vertices 51 . At small scales, the diameter of the polygons is determined by the size of the stress 
shadow around the fractures, which is proportional to the depth to which the fractures 
propagate 51 . The depth of fracturing for small contraction-crack polygons on Earth is dictated by 
the strain rate, the variation of stress with depth, and the rheology of the material in which the 
fractures form 51 . For the PKT, the size of the polygon was likely determined instead by the 
diameter of the tensile stress belt at the surface surrounding the thermal anomaly. The 
propagation of the fractures or rifts into the interior of the region may have been prevented by 
the compressional stresses in the upper lithosphere above the thermal anomaly (Extended Data 
Figs 6-7), which may have had an effect similar to the stress shadows around fractures in small- 
scale polygons. The distribution of KREEP-rich material in the subsurface is poorly constrained, 
while the distribution on the surface is strongly affected by the ejecta of the Imbrium basin and 
the distribution of KREEP-rich maria (the latter of which were controlled in part by the pattern 
of the PKT border structures). An alternative explanation for the pattern of border structures that 
warrants further consideration is that the distribution of KREEP-rich material in the subsurface 
follows a quasi-rectangular pattern. 

Parallels between the PKT and the south polar terrain of Enceladus 

The overall planform of the PKT border structures bears a strong resemblance to that of 
the border structures surrounding the south polar terrain (SPT) of Saturn’s icy moon Enceladus, 
which are also quasi-rectangular in outline 23 . However, as discussed in the main text and 
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expanded upon below, substantial differences exist between these provinces and their in ferred 
evolutionary histories. We emphasize that we do not suggest that the specific processes and 
evolutionary paths of these two regions were identical. Rather, the gross similarities between 
these two provinces on different bodies suggests broad parallels in the processes governing their 
evolution. Here we summarize the basic properties of each province, and then discuss the SPT in 
more detail. 

The PKT on the Moon is a broad area of enhanced surface heat flow as a result of the 
high concentrations of heat-producing elements within the KREEP-rich material 3 ' 10 (Extended 
Data Fig. 8). This compositional anomaly is likely a result of the concentration beneath the 
nearside of the late-stage crystallization products of the lunar magma ocean 34 , including dense 
ilmenite-rich cumulates and KREEP-rich material with high concentrations of U, Th, and K. The 
PKT was the most volcanically active region on the Moon and contains the majority of the mare 
basalt provinces 10 . GRAIL gravity data has revealed the PKT to be surrounded by a quasi- 
rectangular set of magmatic-tectonic structures with straight sides and angular intersections. The 
border anomalies along the northern (Mare Frigoris) and eastern edges of the PKT occur beneath 
maria that are confined within elongated topographic depressions, whereas the border anomalies 
on the western and southern edges of the PKT he adjacent and interior to the topographic step up 
to the highlands. The PKT is characterized by low topography that is largely isostatically 
compensated over long wavelengths. This compensated depression can be explained by a crust 
that is thinner 2 or denser than the surrounding crust, by the presence of denser materials at depth, 
or by a combination of these effects. Thermal annealing of the pore space beneath the PKT due 
to its high heat flow may have increased the bulk density of the crust at depth, which may 
contribute to the low topography 45 . Deeper density anomalies could result from either the 


850 

851 

852 

853 

854 

855 

856 

857 

858 

859 

860 

861 

862 

863 

864 

865 

866 

867 

868 

869 

870 

871 


intrusion of KREEP-rich magma into the lower crust, or from the presence of a remnant of the 
ilmenite-rich cumulates in the upper mantle that may not have fully mixed into the deeper 
interior 37 . Although a thinner crust likely explains most of the observed topography, 
contributions from reduced crustal porosity and the presence of dense materials within or below 
the crust appear likely. 

The SPT on Enceladus is an area of strongly enhanced surface heat flow 26 ' 52 ' 53 (Extended 
Data Fig. 8) as a result of either localized tidal heating or the localized release of global tidal 
heating. The source of this thermal anomaly is thought to be related to the presence of a regional 
liquid water sea or the regional thickening of a global ocean beneath the SPT, which may be a 
result of locally enhanced tidal heating and would itself contribute to the enhanced tidal 
heating 27 ' 54 ' 55 . The SPT is cyrovolcanically active, as revealed by the plume of water vapor and 
icy particles emanating from the parallel “tiger stripe” fractures in the centre of the SPT 24 . 
Cassini Imaging Science Subsystem (ISS) images reveal the SPT to be bound by a quasi- 
rectangular set of tectonic structures with straight sides and angular intersections 24 . These border 
structures occur near the edges of the topographic depression containing the SPT 24 , located either 
at or just above the topographic step leading from the SPT up to the surrounding surface 25 (Fig. 
3c). The SPT is characterized by low topography 2 '^' 56 that is largely isostatically compensated at 
long wavelengths 57 , which is best explained by the presence of a dense subsurface ocean 57 " 58 . 
Depressions in other areas of Enceladus 25 have been explained as a result of the thermal 
annealing of the pore space due to the presence of local thermal anomalies beneath these regions 
in the past 59 . Some contribution to the SPT depression from a reduction of the pore space seems 
probable given the high observed heat flow. However, the large apparent depth of compensation 
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of the SPT indicated by the long-wavelength gravity and topography suggests that the effect of a 
deeper ocean dominates 57 . 

Although there are notable large-scale morphological and geodynamic similarities 
between the PKT and the SPT, there are many differences between these provinces as well. The 
thermal anomaly in the SPT is a result of tidal, rather than radiogenic, heating. Multiple 
mechanisms have been proposed to explain the high heat flux in the SPT, including viscous 
heating in the ice shell 55 and shear heating along fractures 60 . Each of these mechanisms 
ultimately relies on tidal energy from the gravitational interaction of Enceladus with Saturn. 
However, the expected steady-state rate of tidal heating for the present-day eccentricity 61 is not 
sufficient to maintain the observed heat flux within the SPT 26 ’ 53 or the inferred subsurface ocean 
beneath the region 27 . Recent results have revised the lower bounds on the heat flow downward 52 , 
but values remain above the expected steady-state tidal heating unless the dissipation within 
Saturn is higher than expected from theoretical considerations 62 . This discrepancy may be 
explained if the SPT today is in a transient state of high heatflow following an earlier high 
eccentricity phase during which the ocean formed 27 . This scenario implies a time-variable heat 
flux in which the SPT may be cooling today. 

The lack of craters within the SPT 24 suggests an earlier episode of volcanic resurfacing, 
lithosphere recycling 28 ’ 63 , or viscous relaxation of the craters 64 . Each of these scenarios could 
have resulted in a regional thermal anomaly, followed by a period of cooling and contraction of 
the ice throughout the SPT. Globally, substantial lateral and temporal variations in the heat flux 
have been inferred on the basis of high local heat fluxes indicated by the relaxation of craters 64 
and the flexural support of topography 65 . Structures similar in scale and morphology to the SPT 
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on the leading and trailing hemispheres suggest similar activity at those locations in the past 66 , 
further supporting spatial and temporal variability in the thermal state of Enceladus’ ice shell. 

The SPT border structures are each composed of a belt of closely spaced parallel ridges, 
surrounded by an inward (southward) facing scarp 24 ' 67 . The ridge belts likely formed by 
compression 24 ' 25 , though extensional deformation 68 or more complicated scenarios 69 may have 
played a role in the formation of the south-facing scarps. For the compressional interpretation, it 
has been proposed that the tectonics in the SPT was driven by regional thermal expansion 70 , 
which is similar in nature but opposite in sign to what is proposed here for the PKT. At some 
intersections, the border scarps are continuous with fracture belts extending northward from 
120°-angle triple junctions 67 , consistent with an extensional origin for the outer scarp. However, 
the folded terrains confined within the angular corners are indicative of compressional 
deformation 24 . Compressional folding is also observed in the interior of the SPT away from the 
border structures 71 , whereas tensile opening of the “tiger stripe” fractures is required to explain 
the observed volcanic venting 72 . Thus, both compressional and extensional tectonics have been 
active along the border structures and within the interior of the SPT. 

Some structures observed within the SPT are to first order consistent with our model 
predictions for the PKT. The models predict the upper lithosphere within a cooling lithospheric 
cap to be in a state of compression due to its coupling with the contracting lower lithosphere, 
whereas the cap would be surrounded by a belt in which extensional stresses pervade the entire 
lithosphere (Figure 4, Extended Data Figs 6-7). This stress pattern predicts broad compressional 
deformation of the upper lithosphere within the SPT and lithosphere- sc ale normal faulting at the 
margins. However, we emphasize that simple thermal expansion and contraction alone cannot 
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explain the extensive tectonism within the SPT. The extensive tectonic modification and 
resulting large strains may indicate an earlier period of mobile-lithosphere tectonics 28 ' 63 . 

Enceladus is much smaller than the Moon (radii of 252 and 1738 km, respectively). 
Although the SPT is much smaller than the PKT in physical size (-300 km versus -2000 km), 
they are similar in angular size (Fig. 3). Thus, the geometric arguments for the formation of the 
quasi-rectangular PKT border structures due to the intersection of tectonic structures at 120° 
angles on a spherical surface may have relevance to the SPT as well. 

Thus, we suggest that broadly similar geodynamic processes may have been at work in 
the PKT and the SPT. Both regions are characterized by strong thermal anomalies, enhanced 
volcanic activity, and low topography. The quasi-rectangular structures surrounding both 
provinces are consistent with the expected shapes of sets of tectonic structures intersecting at 
120°-angle triple junctions. However, the specific evolutionary paths of the provinces were 
likely substantially different. The source of heat, temporal variations in heat flux, and rheology 
of the lithosphere would have been substantially different in each case. Our current 
understanding of the fonnation and evolution of both structures is incomplete. Nevertheless, both 
provinces highlight the important effect that regional thermal anomalies can have on the regional 
volcanic and tectonic evolution of diverse planetary bodies. 
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1033 Extended Data Table 1. Extension and strain across two border anomalies 


Filter 

p M (kg/m 3 ) 

Anomaly 1 

Anomaly 2 

extension 

strain 

extension 

strain 

isostatic 

3220 

15 km 

0.11 

13 km 

0.08 

equal amplitude 

3220 

12 km 

0.09 

10 km 

0.07 

isostatic 

3500 

1 1 km 

0.09 

10 km 

0.06 

equal amplitude 

3500 

1 1 km 

0.08 

9 km 

0.06 

mare-crust only 

3220 

10 km 

0.07 

8 km 

0.05 

crust-mantle only 

3220 

16 km 

0.12 

18 km 

0.12 
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Extended Data Figure 1. Comparison of the GRAIL gravity gradients with proposed 
Procellarum basin ring structures, a, Bouguer gravity gradients (Eotvos) in a Lambert 
azimuthal equal-area projection of the nearside of the Moon, b, Muted gravity gradients overlaid 
with mapped mare shorelines and scarps (dots) and wrinkle ridges (lines) (modified from Figure 
1 of Whitaker 5 ). 
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Extended Data Figure 2. Filters applied during the crustal thickness modelling. Filters 
were applied during the calculation of the relief along the crust-mantle interface (solid) and 
mare-crust interface (dashed) for cases in which the relief along the two interfaces was either a, 
isostatic prior to mare loading or b, equal and opposite in amplitude. The filter in b imposes the 
isostatic condition from degrees 1 to 3, with a linear transition to the equal amplitude filter from 
degrees 3 to 10. Both filters apply a cosine taper from degrees 125 to 150. The mare-crust filter 
is shown for illustration purposes only. In practice, the relief along the mare-crust interface was 
calculated using the residual Bouguer anomaly after the calculation of the crust-mantle interface 
relief (equivalent to using the filter shown with the original Bouguer gravity). 
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Extended Data Figure 3. (next page) Predicted thicknesses of the crust and maria and 
average cross-sections across two of the anomalies. Predicted thickness of the maria (left) and 
underlying feldspathic crust (middle), and cross-sections of the modeled structure of anomaly 1 
(right, top) and anomaly 2 (right, bottom) showing the variation in the thickness of the mare 
(dark gray) and feldspathic crust (light gray). Models are for the cases of a-d, isostatic relief 
along the two interfaces prior to mare infilling with a mantle density of 3220 kg/m 3 , e-h, equal- 
amplitude relief along the two interfaces with a mantle density of 3220 kg/m 3 , i-1, isostatic relief 
along the two interfaces prior to mare infilling with a mantle density of 3500 kg/m 3 , m-p, equal- 
amplitude relief along the two interfaces with a mantle density of 3500 kg/m 3 , q-t, all gravity 
anomalies at degrees >10 ascribed to relief on the mare-crust interface, and u-x, all gravity 
anomalies at degrees >10 ascribed to relief on the crust-mantle interface. 
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Extended Data Figure 4. Temperature evolution at a depth of 25 km. Temperatures are 
shown for a, KREEP-rich material concentrated at the base of the crust and b, KREEP-rich 
material distributed throughout the crust. The temperature at the base of the crust outside the 
PKT is shown for comparison. The period between 4.0 and 3.0 Ga that is the focus of the stress 
modeling is indicated by the shaded box. 
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Extended Data Figure 5. Predicted changes in temperature relative to areas outside the 


PKT and absolute temperature change between 4.0 and 3.0 Ga. Results are shown for cases 


with a,b, KREEP concentrated below the crust and c,d, KREEP distributed throughout the crust. 


The PKT is centred on the pole at the left side of the panels. The region shown in Extended Data 


Figs 6-7 (encompassing 90° of arc extending radially outward from the centre of the PKT and 


downward to a depth of 50 km) is outlined in black. 
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Extended Data Figure 6. Predicted lithospheric stresses and magma ascent for the case of 
10 km of KREEP at the base of the crust. Cross-sections show a, the in-plane horizontal 
stresses (radial to the centre of the PKT, the far-field stress profile was subtracted to calculate the 
relative stress), b, the difference between the in-plane horizontal stress and the vertical stress, c, 
the magma ascent criteria, and d, the deviatoric stress. The magma ascent criteria reveals 
portions of the crust in which the horizontal stresses are tensile relative to the vertical stresses to 
permit the formation of vertical dikes (dark gray), where the vertical stress gradient is more 
favorable to magma ascent than the lithosphere far from the PKT (light grey), where magma will 
rise unassisted by other factors such as pressurized magma chambers (red), and where none of 
the criteria are satisfied (diagonal lines). 
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Extended Data Figure 7. Predicted lithospheric stresses and magma ascent for the case of 


10 km of KREEP basalt distributed uniformly through a 40-km-thick crust. All panels are 
as in Extended Data Fig. 6. 
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Extended Data Figure 8. Additional comparisons of Procellarum KREEP terrane to the 
Enceladus south polar terrain, a, The PKT is characterized high heat flow as a result of the 
enhanced abundances of radioactive elements 3 (represented by the concentration of thorium 4 ), b, 
The border structures of the SPT as revealed by Cassini ISS images 24 also trace a quasi- 
rectangular pattern enclosing a region of c, elevated brightness temperatures and enhanced heat 
flow 26 . All maps are in a simple polar projection. In all panels, the circle corresponds to an 
angular diameter of 180° of surface arc, divided into 10° increments. 


